A functional of the one-body-reduced density matrix derived from the homogeneous 

electron gas: Performance for finite systems 
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An approximation for the exchange-correlation energy of reduced-density-matrix-functional the- 
ory was recently derived from a study of the homogeneous electron gas (N.N. Lathiotakis, N. Helbig, 
E.K.U. Gross, Phys. Rev. B 75, 195120 (2007)). In the present work, we show how this approx- 
imation can be extended appropriately to finite systems, where the Wigner Seitz radius r s , the 
parameter characterizing the constant density of the electron gas, needs to be replaced. We apply 
the functional to a variety of molecules at their equilibrium geometry, and also discuss its perfor- 
mance at the dissociation limit. We demonstrate that, although originally derived from the uniform 
gas, the approximation performs remarkably well for finite systems. 
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I. INTRODUCTION 



Density-functional theory (DFT) is one of the most 
widely applied tools for electronic structure calculations. 
While it successfully describes systems ranging from 
atoms and molecules to solids, present day DFT approxi- 
mations fail to describe a class of systems generally called 
strongly correlated. For these systems, recent calcula- 
tions have fueled the hope that reduced-density-matrix- 
functional theory (RDMFT) can cure the problem. 1-3 
Within RDMFT, the one-body reduced density matrix 
(1-RDM), 7(r,r'), is used as the basic variable in anal- 
ogy to DFT, where that role is reserved for the electronic 
density. The theorem of Gilbert , which is an extension 
of the Hohenberg-Kohn theorem to non-local external 
potentials, guarantees that the ground-state expectation 
value of any observable of a quantum mechanical system 
is a unique functional of the ground-state 1-RDM. In par- 
ticular, the total energy E tot of a system of N electrons 
moving in an external local potential V(r) can be written 
in terms of the ground-state 7, as 



Etotil] = // d 3 rdV5(r-r') 



-Iv r 2 + U(r) 



7 (r, r') 



+ \Jj M r ^^K E M (1) 

(atomic units are used throughout). The first term con- 
tains the kinetic and external energies and is a sim- 
ple functional of 7. The fact that the kinetic energy 
can be written explicitly in terms of the ground-state 1- 
RDM is a great advantage of RDMFT compared to DFT. 
The energy contribution associated with the electron- 
electron interaction can be cast into two terms, the di- 
rect Coulomb energy (or Hartree) term which is again an 
explicit functional of 7, and the remaining contribution 



which is called exchange-correlation (xc) energy. The 
xc energy is an unknown functional of the 1-RDM and 
needs to be approximated in practice. Contrary to DFT, 
however, it does not contain any kinetic energy contri- 
butions but is solely given as the difference between the 
full Coulomb interaction and the Hartree energy. Sev- 
eral approximations for the xc energy have been intro- 
duced so far 1-3,5-16 , the great majority of them being 
implicit functionals of the 1-RDM. They depend explic- 
itly on the natural orbitals ipj and the corresponding oc- 
cupation numbers rij, i.e. the eigenfunctions and the 
eigenvalues of 7, which are given by 



/ 



dV 7 (r,rVi(r') = ?w(r). 



(2) 



Applications of different RDMFT functionals for the cal- 
culation of the dissociation of molecules 2,17 the ionization 
potential 18-21 , or the fundamental gap 21 ' 22 , have been re- 
ported. Most approximate RDMFT functionals can be 
written in the form 

£ xc [7] = £*c [{rij},{<pj}] = 



3,1=1 



(3) 



with some function f(nj,ni) which distinguishes these 
approximations from the Hartree-Fock approximation. 

The first approximation of the kind in Eq. (3) was 
introduced by Miiller 5 using the function f(rij,rii) = 
^Jnjni. Buijse and Baerends rederived the same approx- 
imation from modeling exchange and correlation holes. 1 
The Miiller functional overcorrelates in all systems it was 
applied to. 23-28 Interestingly, for diatomic molecules, it 
reproduces the dissociation into two independent atomic 
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fragments correctly. Goedecker and Umrigar 6 (GU) con- 
sidered a functional with the same function / but exclud- 
ing explicitly the self-interaction terms, j = I, in Eq. (3) 
and in the direct Coulomb term in Eq. (1). They used 
their functional in a minimization procedure and deter- 
mined the natural orbitals and the corresponding occu- 
pation numbers. In this way, they found correlation en- 
ergies, for small atomic and molecular systems, that are 
in very good agreement with the exact ones. However, it 
was found later that the GU functional fails to describe 
the dissociation of diatomic molecules correctly 23,24 . In 
an attempt to correct the overcorrelation of the Muller 
functional, while keeping the good description in the dis- 
sociation limit, Gritscnko et al. 2 proposed a hierarchy 
of three levels of repulsive corrections. All three correc- 
tions distinguish between weakly and strongly occupied 
orbitals, the former being orbitals with occupation num- 
bers close to zero, the latter having occupations close to 
one. The resulting functionals are called BBC1, BBC2, 
and BBC3. For the first two approximations the function 
/ is given by 

J-y/fTj m , j ^ I, and 
j, I weakly occupied, 
y'rij ni , otherwise, 

(4) 

{-^njnj, j ^ I, and 
j, I weakly occupied, 
Hj m, j ^ I, and 
j, I strongly occupied, 
^Jn.j ni , otherwise. 

(5) 

In the BBC3 functional, the anti-bonding orbital is 
treated as strongly occupied. Additionally, the self- 
interaction terms are removed as in the GU functional, 
except for the pair of bonding and anti-bonding orbitals. 
Gritsenko et al used the BBC functionals to calculate 
the dissociation curves of diatomic molecules. They con- 
cluded that the BBC3 functional is very accurate in the 
description of these systems both at the equilibrium ge- 
ometry and the dissociation limit. Two other function- 
als, derived from a cumulant expansion of the second or- 
der density matrix, with a final form that is very similar 
to the BBC functionals were introduced by Piris. 14 We 
refer to these functionals as Piris natural orbital func- 
tionals. The first approximation, PNOFO, is identical 
to the BBC1 functional apart from the self- interaction 
terms which are removed in the same way as in the GU 
approximation. In the second functional, PNOF, an ad- 
ditional term is included to avoid occupation numbers 
which arc identical to zero or one. PNOFO and PNOF 
coincide for two electron systems. The BBC as well as 
PNOF and PNOFO functionals were evaluated recently 
for a large set of molecular systems and were proven to 
be quite accurate in reproducing the correlation and the 
atomization energies of these systems. 25 

For the application of 1-RDM functionals to periodic 
systems, the homogeneous electron gas (HEG) is an im- 



portant prototype system. Also, as far as size is con- 
cerned, the HEG and small atomic and molecular systems 
are two opposite extremes. For the HEG, the GU and 
Muller functionals are identical since the self-interaction 
terms vanish. Similarly, all terms that include a spe- 
cial treatment of single orbitals vanish. As a result, the 
BBC3 functional coincides with BBC2, and PNOFO with 
BBC1 in this special case. Cioslowski and Pernal 26 as 
well as Csanyi and Arias 27 applied the Muller functional 
to the HEG. As for finite systems, the correlation energy 
is overestimated, actually by a factor of about two in 
the high density regime. In the low density region, the 
Muller functional fails completely to reproduce the limit 
of zero correlation. 26-28 Csanyi and Arias also introduced 
a different functional starting from a tensor expansion of 
the second-order matrix. Unfortunately, it fails to de- 
scribe the electronic correlation of the HEG in both the 
dense and dilute limits. A more successful functional for 
the HEG is the one proposed by Csanyi, Goedecker and 
Arias (CCA) 7 , which reproduces relatively accurately the 
correlation energy for the dense HEG. Recently, we ap- 
plied the BBC1 and BBC2 functionals to the HEG 28 . 
We showed that these functionals offer a better descrip- 
tion of the correlation of the HEG over the whole range 
of densities than any other of the discussed functionals. 
Both the BBC1 and the BBC2 functionals overcorrelate 
slightly for high densities and undercorrelate for low den- 
sities. The crossover is at around r s = 0.5 where these 
functionals perform best. Additionally, they produce a 
finite discontinuity in the momentum distribution at the 
Fermi energy, resembling a feature of the exact theory. 
Unfortunately, the size and the dependence on the den- 
sity of this discontinuity are not in agreement with the 
exact result. 



II. THE BBC++ FUNCTIONAL 

In an attempt to improve over the BBC functionals 
for the HEG, we introduced a modification to the BBC1 
functional 28 . It consists of the introduction of a param- 
eter into the function / fitted for each value of r s to 
reproduce the exact correlation energy of the HEG. Two 
choices for this parameter were made. We concluded that 
a reasonable way to generalize BBC1 is to introduce a 
function s(r s ) multiplying the function / when both or- 
bitals are weakly occupied. The corresponding function 
/ then reads 

!—s(r s ) y/rij ni , j^/,and 
j, I weakly occupied, 
y/rij ni , otherwise. 

(6) 

We call this functional s-functional. Note that, for the 
HEG, if one assumes plane-wave natural orbitals, no spe- 
cial care is necessary for the j = I terms since these terms 
vanish. The s-functional reproduces, by construction, the 
exact correlation energy of the HEG. It also improves the 
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consequently a, are only known at the solution point, i.e. 
when the optimal 7 is known. Thus, s(a) has to be de- 
termined sclf-consistently during the minimization proce- 
dure. Starting from a trial value for s one minimizes the 
energy with respect to 7, calculates a and feeds the cor- 
responding value for s back into the functional. As the 
BBC++ functional coincides with the s-functional for 
the HEG, the self-consistent determination of s(a) has 
to yield the correct s(r s ) in this case. Using the imple- 
mentation for the HEG presented in Ref. 28, we verified 
that this is indeed the case. We also found that the self- 
consistent determination of s converges for all the finite 
systems we studied. 
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FIG. 1: Dependence of the fitting parameter s of the s- 
functional 28 on the ratio a of the correlation over the kinetic 
energy for the HEG. Two different functions s are plotted 
for two sets of diffusion Monte-Carlo data which s was fit- 
ted to reproduce: From Ceperley and Alder 29 and Ortiz and 
Balone 30 . 



calculated momentum distribution compared to BBC1 
and BBC2. However, the momentum distribution still 
deviates from the exact one 28 . 

Being derived from the study of the HEG, the s- 
functional is directly applicable to metallic systems, 
where a value of r s can be associated. The application 
to other systems is not straightforward since one has to 
overcome the r s dependence. One possibility, in the spirit 
of the local density approximation of DFT, is to express s 
as a function of the local electronic density. Hence, s be- 
comes a function of the space coordinate r and therefore 
should be properly included in the integrals in Eq. (3). 
A reasonable choice that preserves the symmetry of the 
integrations over r and r' is to multiply the integrand in 
Eq. (3) by y/ s(n(r)) s(n(r')) or other possible averages 
and keep the BBC1 form for the function /. Work inves- 
tigating the performance of the resulting functional is in 
progress. 

In the present article, we propose an alternative way to 
circumvent the r s dependence. We refer to the resulting 
functional as the BBC++ functional. In this approxima- 
tion, the xc terms retain the simple form of the exchange 
integrals over Gaussian or Slater type orbitals. The idea 
is to establish the dependence of s on a quantity a which, 
contrary to r s , is meaningful for all systems, finite and 
periodic, metallic and insulating. We choose the ratio 
of the correlation energy over the Hartrcc-Fock kinetic 
energy for this quantity a. For the HEG, a depends on 
the density parameter r s and, since the function a(r s ) 
is strictly monotonic and therefore invertible, the depen- 
dence s(a) can be established. The resulting function 
s(a) for the HEG is shown in Fig. 1. The function s(r s ), 
which is necessary in the calculation of s(a), is given in 
Ref. 28. The difficulty for the application of the BBC++ 
functional lies in the fact that the correlation energy, and 



III. RESULTS FOR FINITE SYSTEMS 

For the application of 1-RDM functionals to finite sys- 
tems, we used the implementation introduced in Ref. 12 
which relies on the GAMESS computer program 31 for the 
calculation of the one- and two- electron integrals. In the 
present work, we employed the cc-pVTZ basis set 32 for 
all systems apart from the He atom for which the cc- 
pVQZ set 33 was used. Depending on the system, these 
basis sets contain 30 to 50 basis functions. We always 
made full use of the size of the basis set, optimizing as 
many natural orbitals as there were basis functions avail- 
able. For the BBC3 functional, we used the form that 
respects the possible degeneracies of the bonding and 
anti-bonding orbitals. 25 The total energies resulting from 
this full minimization with respect to the natural orbitals 
and occupation numbers for several atoms and diatomic 
molecules are given in Table I. We compare all our results 
to total energies obtained from a coupled-cluster-singles- 
doubles-triples (CCSD(T)) 34 calculation using the Gaus- 
sian 03 computer program 35 with the same basis sets as 
used in the RDMFT calculation. For the systems consid- 
ered here, the BBC++ functional yields slightly better 
total energies than GU but docs not reach the accuracy 
of the BBC3 and Piris functionals. However, it performs 
significantly better than the BBC1 approximation that 
it was derived from, except for the H2 molecule. The re- 
pulsive correction to the BBC1 functional increases with 
increasing number of electrons in the system. Overall, for 
small finite systems, the BBC++ functional performs re- 
markably well considering that it was originally tuned to 
be exact for the HEG. 

In Fig. 2, we plot the dissociation curve for the H2 
molecule. As already mentioned, among all the sys- 
tems we studied, H2 is the only case where the BBC++ 
functional does not improve over BBC1. Nevertheless, 
BBCH — h reproduces the dissociation curve surprisingly 
well giving a qualitatively correct optimal 1-RDM with 
four occupations equal to 0.5 (per spin) in the limit of 
large distance. The slight decrease in the total energy 
at around 2.5A separation, which can be seen in Fig. 2, 
is a pathology of the BBC++ functional. It originates 
from the dramatic increase in static correlation leading 
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Functional 


He 


Be 


H 2 


LiH 


Nc 


HF 


A 


Hartree-Fock 


2.8615 


14.5729 


1.1330 


7.9868 


128.5320 


100.0589 


0.119 


Miiller 


2.9143 


14.7471 


1.1905 


8.1167 


128.9168 


100.5185 


0.090 


GU 


2.8980 


14.6456 


1.1660 


8.0478 


128.8260 


100.3692 


0.019 


BBC1 


2.9042 


14.6692 


1.1746 


8.0642 


128.8471 


100.4105 


0.035 


BBC3 


2.9005 


14.6490 


1.1705 


8.0444 


128.8014 


100.3373 


0.011 


PNOF 


2.8925 


14.6234 


1.1593 


8.0310 


128.7876 


100.3166 


0.013 


PNOFO 


2.8925 


14.6236 


1.1593 


8.0309 


128.7840 


100.3133 


0.014 


BBC++ 


2.9035 


14.6574 


1.1768 


8.0572 


128.7862 


100.3526 


0.018 


CCSD(T) 


2.9025 


14.6186 


1.1724 


8.0227 


128.8054 


100.3401 


0.000 


s 


1.202 


2.108 


0.573 


1.674 


3.158 


2.823 




a 


-0.0147 


-0.00580 


-0.0388 


-0.00882 


-0.00198 


-0.00344 





TABLE I: Absolute total energies (atomic units) of atoms and diatomic molecules calculated with different 1-RDM functionals 
and the average absolute deviation A from the reference energies obtained with the CCSD(T) 34 method. For the latter, we 
used the Gaussian 03 program 35 for the same basis sets. In the last two rows, the values of the optimal values of s for the 
BBC+- 1- functional and the corresponding values of a are also given. 
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FIG. 2: Dissociation curve of the H2 molecule with different 
functionals of the 1-RDM. 

to a large negative value for a. While a has a value of 
-0.0388 at the equilibrium distance, it becomes -0.335 at 
a distance of 4 A. As we can see from Fig. 1, this results 
in a negative value for s driving the BBC++ towards 
the Miiller functional. Consequently, the functional over- 
correlates and leads to a decrease in the total energy 
with increasing distance for the H2 molecule. It is worth 
mentioning that the dissociation of the H2 molecule is a 
rather difficult case for DFT functionals 36 ' 37 . 

Finally, a second pathology of the BBC++ functional 
is its obvious size inconsistency: Consider a system con- 
sisting of two independent sub-systems, for example two 
finite systems at a large distance. If the two sub-systems 
are identical, the a of the total system is equal to the 
values of each sub-system. In the extreme case, how- 
ever, where one of the sub-systems is much larger than 
the other, the common value for a is completely deter- 
mined by the larger sub-system. That is, since the two 
systems arc independent, the functional, when applied 
to the composite system, gives a different result for the 



smaller sub-system than when applied to the sub-systems 
independently. 



IV. CONCLUSION 

We presented a RDMFT functional, which we call 
BBC++, based on an idea to circumvent the dependence 
on the density parameter r s of functionals derived from 
the homogeneous electron gas. This idea is applied to the 
s-functional introduced in Ref. 28. We apply BBC++ in 
the calculation of correlation energies of small atomic and 
molecular systems and show that its performance is sat- 
isfactory. We also discuss pathologies of this functional, 
with the most important being its size inconsistency. 

Despite its pathologies, the BBC++ functional repre- 
sents an important step in the development of 1-RDM 
functionals. It is a successful attempt to apply approx- 
imations originally developed for the HEG to finite sys- 
tems. In other words, within RDMFT it is possible to de- 
velop functionals that perform equally well for extended 
systems, like the HEG, as well as small atomic and molec- 
ular systems. The present work serves as an initiative for 
the development of better approximations based on the 
HEG and, furthermore, their application to finite systems 
in the future. 
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